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Annual  Summary  Report 
“Quantum  Monte  Carlo  for  Molecules” 

Principal  Investigators: 

William  A.  Lester,  Jr. 

Peter  J.  Reynolds 

Materials  and  Molecular  Research  Division 
Lawrence  Berkeley  Laboratory 
University  of  California 
Berkeley,  California  94720 

Description  of  Problem  and  Approach 

Introduction. 

Many-body  problems  in  physics  are  often  treated  by  a  Monte  Carlo 
approach  [l].  The  Monte  Carlo  method  is  statistical  in  nature.  Based  on  the  gen¬ 
eration  of  "random"  numbers  or  "coin  tosses,"  it  derives  its  name  from  a  city 
famous  for  the  random  numbers  embodied  in  its  games  of  chance.  Thus  it 
perhaps  easy  to  imagine  using  the  Monte  Carlo  method  for  treating  inherently 
statistical  models,  or  even  for  numerical  integration  [2]:  it  is,  however,  less  obvi¬ 
ous  how  to  solve  many-body  problems  with  it.  Nevertheless,  many  such  prob¬ 
lems  are  readily  treated  by  Monte  Carlo. 

Of  particular  interest  to  molecular  physics  are  the  quantum  mechanical 
Monte  Carlo  (or  QMC)  methods^.  What  we  mean  by  QMC  is  a  Monte  Carlo 
procedure  for  solving  the  Schrodinger  equation  statistically  by  the  simulation  of 
an  appropriate  random  process.  The  formal  similarity  of  the  Schrodinger  equa¬ 
tion  with  a  diffusion  equation  allows  one  to  calculate  quantum  mechanical  expec¬ 
tation  values  as  Monte  Carlo  averages  over  an  ensemble  of  random  walks.  The 


2 


method  is  procedurally  quite  simple.  As  a  result,  QMC  provides  an  attractive 
alternative  to  the  conventional  variational  and  perturbation-theoretic  techniques 
used  in  physics  and  chemistry.  — '  ,■  / ' u 


Past  Work 

We  have  applied  QMC  sucessfully  to  the  calculation  of  the  total  energy  of  a 
number  of  molecular  systems  [4-7].  In  every  case  we  have  achieved  very  high 
accuracy  compared  with  experimentally  inferred  and  exact  results  (where  avail¬ 
able),  as  well  as  ab  initio  configuration  interaction  calculations.  In  most  cases, 
90-100%  of  the  correlation  energy  has  been  obtained. 

In  the  past  few  years  our  research  program  has  begun  to  make  a  more 
significant  contribution  [5-7]  by  obtaining  quantities  of  more  physical  interest 
than  total  energies.  For  example,  much  of  chemistry  takes  place  predominantly 
in  the  valence  electrons  of  a  system.  Thus  the  quantities  of  interest  are  usually 
small  differences  of  large  total  energies.  Examples  of  such  differences  include 
binding  energies,  electron  affinities,  reaction  barriers,  and  level  splittings.  We 
have  performed  calculations  for  each  of  the  above-mentioned  quantities.  Such  cal¬ 
culations  are  a  far  more  difficult  task  for  Monte  Carlo,  since  a  small  statistical 
uncertainty  (e.g.  of  as  little  as  0.1°o)  in  the  separate  total  energies  can  mask  the 
sought-after  energy  difference.  To  reduce  the  statistical  error  to  the  level  needed 
by  “brute  force”  is  costly  in  computer  time,  as  the  standard  deviation  decreases 
only  as  ( CPU  time  )~* .  Algorithmic  developments,  such  as  differential  QMC  [8] 
hold  promise  for  reductions  in  variance  through  correlated  sampling  techniques. 
Another  approach  is  based  on  the  variance  reduction  achieved  as  a  trial  wave 


function  'I'  j  approaches  the  true  eigenfunction.  To  take  advantage  of  this  latter 
approach  we  have  developed  an  iterative  procedure  for  improving  T  [7].  We 
are  currently  working  on  an  improved  iterative  procedure,  as  well  as  on  different 
forms  of  correlation  functions,  and  on  methods  of  Monte  Carlo  optimization  of 
the  correlation  parameters. 

Using  QMC,  we  have  recently  studied  points  along  the  reaction  coordinate  of 
the  H  +  H2  exchange  reaction  [6] .  Particular  emphasis  has  been  placed  on  the 
saddle-point  geometry,  for  which  Liu  [0]  has  performed  the  most  extensive  Cl  cal¬ 
culation  to  date.  The  bound  for  the  barrier  height  which  we  obtained  by  QMC  is 
0.16  kcal/mole  below  Liu’s  bound,  and  probably  lies  within  0.1  kcal/mole  of  the 
exact  answer.  In  addition,  we  were  able  to  obtain  these  results  with  only  single¬ 
determinant  trial  functions,  and  a  basis  set  expansion  at  only  the  double-zeta 
level.  The  nodes  (see  below),  which  are  important  in  determining  the  correct 
energy,  prove  to  be  quite  insensitive  to  basis  set  beyond  the  double-zeta  level. 

Until  recently,  QMC  applications  have  been  limited  to  calculations  of  ener¬ 
gies  of  ground  states  and  lowest  states  of  a  given  symmetry.  As  an  example  of 
the  latter,  we  calculated  the  singlet-triplet  splitting  in  methylene  [7],  The  two 
states  studied  are  both  lowest-energy  states  of  their  respective  symmetries.  In 
these  studies,  the  accuracies  obtained  with  QMC  have  been  comparable  to  the 
best  achieved  by  conventional  ab  initto  methods.  In  section  II,  where  we  discuss 
progress  during  the  current  year,  we  include  our  recent  effort  in  extending  QMC 
to  excited  states  of  the  same  symmetry  as  the  ground  state.  In  that  section,  we 
also  report  on  calculations  of  properties  unrelated  to  the  energy  and  of  energy 


1  -Y -Y-YsY  V A  -Y-Y-Y 


l  fc-|  H«1 


& 


•VV 


Y  Y  VY",  »jY  "Y'. ..  ^ 


gradient  calculations.  However,  before  going  into  detail,  in  what  follows  we 
describe  the  theoretical  background  of  QMC. 

QMC  Theory 

Briefly,  one  simulates  a  quantum  molecular  system  by  allowing  it  to  evolve 
under  the  time-dependent  Schrodinger  equation  in  imaginary  time.  It  is  easy  to 
show  [4]  that  the  use  of  imaginary  time  causes  the  system  to  approach  a  station¬ 
ary  state  which  is  the  lowest  state  of  a  given  symmetry.  Many  properties  may 
then  be  “measured”  as  averages  over  the  resulting  equilibrium  distribution. 

By  writing  the  imaginary-time  Schrodinger  equation  with  a  shift  in  the  zero 
of  energy  as 

d^dt'1'  [ET-V(R)}*(R,t)  ,  (1) 

we  see  that  it  may  be  interpreted  as  a  generalized  diffusion  equation.  The  first 

term  on  the  right-hand-side  is  the  ordinary'  diffusion  term,  while  the  second  term 
is  a  position-dependent  rate  (or  branching)  term.  For  an  electronic  system. 
D  =fr /2me ,  R  is  the  three-N  dimensional  coordinate  vector  of  the  N  electrons, 
and  V{&)  is  the  Coulomb  potential.  Since  diffusion  is  the  continuum  limit  of  a 
random  walk,  one  may  simulate  Eq.  (1)  with  the  function  'k  (note,  not  'I'2)  as  the 
density  of  “walks”.  The  walks  undergo  an  exponential  birth  and  death  as  given 
by  the  rate  term.  This  connection  between  a  quantum  system  and  a  random 
walk  was  first  noted  by  Metropolis,  who  attributes  it  to  Fermi  [10]. 

The  steady-state  solution  to  Eq.  (1)  is  the  time-independent  Schrodinger 
equation.  Thus  we  have  ^(Z?  .t  )—*<t>{R.  ).  where  O  is  an  energy  eigenstate.  The 


value  of  Et  at  which  the  population  of  walkers  is  asymptotically  constant  gives 
the  energy  eigenvalue.  Early  calculations  employing  Eq.  (1)  in  this  way  were 
done  by  Anderson  on  a  number  of  one-  to  four-electron  systems  [l  1  ]. 

Unfortunately,  in  order  to  treat  systems  larger  than  two  electrons,  the  Fermi 
nature  of  the  electrons  must  be  taken  into  account.  The  antisymmetry  of  the 
eigenfunction  implies  that  must  change  sign;  however,  a  density  (e.g.  of  walk¬ 
ers)  cannot  be  negative.  To  handle  this,  Anderson  made  simplifying  assumptions 
about  the  positions  of  the  nodes.  His  method  was  ad  hoc  ,  and  not  readily  gen- 
eralizable.  Another  method  which  imposes  the  antisymmetry,  and  at  the  same 
time  provides  more  efficient  sampling  (thereby  reducing  the  statistical  “noise”),  is 
importance  sampling  with  an  antisymmetric  trial  function  'I' T  (see  e.g.  Ref.  4). 
The  zeroes  (nodes)  of  'I'j  become  absorbing  boundaries  for  the  diffusion  process: 
this  maintains  the  antisymmetry.  The  additional  boundary  condition  that 
vanish  at  the  nodes  of  ^  T  is  the  fixed-node  approximation  [4,12].  The  magni¬ 
tude  of  the  error  thus  introduced  depends  on  the  quality  of  the  nodes  of  ^j-iB ), 
and  vanishes  as  approaches  the  true  eigenfunction.  To  the  extent  that  'k  T  is 
a  good  approximation  of  the  wave  function,  the  true  eigenfunction  is  almost  cer¬ 
tainly  quite  small  near  the  nodes  of  'k  T .  Thus  one  expects  the  fixed-node  error 
to  be  small  for  reasonable  choices  of  'k  T  .  Work  on  a  number  of  systems  has 
borne  this  out  [4-7,13,14].  In  addition,  this  error  is  variationally  bounded. 

To  implement  importance  sampling  and  the  fixed-node  approximation.  Eq. 
(1)  is  multiplied  by  'P  j- ,  and  rewritten  in  terms  of  the  new  probability  density 
/  (R  ,t  )='i>  r  (B  )V(B  J  )•  The  resultant  equation  for  /  (/?  ,t  )  may  be  written 


-^-  =  Z>V-7  +{ET-EL(R)}f  -DV-lfFQ(R)}  .  (2) 

The  local  energy  EL(R ),  and  the  “quantum  force”  Fq(R)  are  simple  functions 
of  given  by 


El(R)=H*t(R)/*t(R)  , 


Fq(R)=2V*t(R)/*t(R)  .  (3b) 

Equation  (2),  like  Eq.  (1),  is  a  generalized  diffusion  equation,  though  now  with  the 

addition  of  a  drift  term  due  to  the  presence  of  Fq  . 

In  order  to  perform  the  random  walk  implied  by  Eq.  (2)  we  use  a  short-time 
approximation  to  the  Green’s  function.  The  Green’s  function  is  used  to  evolve 
the  distribution  forward  in  time,  i.e.  f  {R  ,t  )—>  f  {R1  ,t +t).  This  process  is 
iterated  to  large  t  .  Such  an  approach  becomes  exact  in  the  limit  of  vanishing 
time-step  size,  r.  Asymptotically,  /  (R  ,t  )—  /00( R  )='!' T  (R  )<t>(R  ). 

The  function  4>{R  )  is  the  lowest-energy  eigenfunction  of  the  Schrodinger 
equation  for  the  imposed  set  of  nodes.  Although  neither  this  function  nor  fx  is 
known  analytically,  we  can  nevertheless  sample  desired  quantities  from  the  equili¬ 
brium  distribution.  Averages  taken  with  respect  to  the  distribution  fx  are 
known  as  mixed  averages.  For  example,  sampling  a  quantity  A  in  equilibrium 
after  A'  samples  gives  the  average  (in  the  limit  of  large  .V) 
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/#  T{R  )<?(/?)  A  dR 

J*r(£)o(E)  dR 

or  in  abbreviated  Dirac  notation  (with  the  normalization  absorbed), 

<A  >  /x=  <VT  |-4  |  q>  .  (5) 

On  the  other  hand,  the  correct  expectation  value  of  A  ,  for  a  state  d.  is 

<0  |  -4  |  0>.  In  computing  the  energy,  or  any  property  for  which  o  is  an  eigen¬ 
state,  there  is  no  difference  between  these  two  averages.  This  follows  since  the 
eigenvalue  can  be  taken  out  of  the  integral  in  the  numerator  of  Eq.  4.  In  particu¬ 
lar,  to  compute  the  energy  one  samples  the  local  energy  EL  ( R_ ).  Then 

fo(R  )♦  T  (&  )[^  tUL  )I{  *  T  UL  )]  dZ 
/o(£)4'r(/?  )  dR 

—  <o  |  H  |  >  — £'0  ,  (6) 

where  E0  is  the  eigenvalue  corresponding  to  the  state  o.  The  last  equality  fol¬ 
lows  upon  noting  that  //  is  Hermitian.  and  thus  can  operate  to  the  left.  This 
ends  our  overview  of  QMC  theory. 

Progress  in  Current  Year 

Molecular  Properties.  For  expectation  values  of  quantities  whose  opera¬ 
tors  do  not  commute  with  // ,  the  mixed  average  of  Eq.  5  is  only  approximate. 
One  suspects  that  the  mixed  average  is  in  some  sense  "half-way"  between  the 
exact  expectation  value  (with  respect  to  o)  and  the  variational  expectation  value, 
taken  with  respect  to  the  trial  wave  function,  i.e.  <  'k  T  |  .4  |  T  > .  Taken 
literally,  this  implies  t  hat  <0  [  .4  |  o>  =  2<  'l'  T  |  .4  j  o>  -  <  4*  T  |  .4  |  'I'  T  > . 
This  result  can  be  formalized  through  the  following  argument.  The  trial  function 
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'I'j-,  if  it  is  good,  differs  from  only  by  a  “small”  function  A,  i.e.  ^=^r+A. 
Then 

<4>\  A  |  4>>  =  <Vt  |  A  |  4»  +  <A  j  A  \  4>> 

^C'l'r  |  .4  |0>  +  <A|A  |^r> 

=2<'I'r  |  A  |  4>>  -  <VT  |  A  |  *T  >  .  (7) 

The  first  equality  follows  on  expanding  the  4>  bra.  The  approximation  in  the  next 

line  occurs  on  expanding  the  <t>  ket  in  the  second  term,  and  dropping  the  resulting 
term  of  order  A2.  Finally,  A  is  re-expanded;  A  is  assumed  an  observable,  and 
hence  Hermitian.  This  gives  an  approximate  formula  for  expectation  values 
taken  solely  with  respect  to  <j>  from  just  mixed  and  variational  averages.  The 
above  argument  ignores  the  different  normalizations  implicit  in  the  different 
terms.  However,  it  is  easy  to  demonstrate  that  Eq.  7  divided  by  <0  |  0>  differs 

from  2<'kr  |.4  |0>/<'kr  |  4»-  <V  T  |  A  \  V  T  >  /  <V  T  |  't'j-  >  by  terms 

of  only  O  (A2).  This  gives  the  desired  result.  It  is,  however,  difficult  to  know 
how  significant  it  is  to  drop  terms  of  order  A2.  Thus,  it  is  desirable  to  be  able  to 
sample  exactly  from  the  distribution  j  <$>  |  2.  This  can  be  done,  though  with  some 
changes  to  the  usual  QMC  algorithm.  The  distribution  J ^  must  be  weighted 
locally  by  0{R)/'l'j{R).  This  quantity  is  essentially  the  asymptotic  number  of 
survivors  of  the  local  configuration  /?  [15].  Thus,  algorithmically,  one  must  fol¬ 
low  each  configuration  into  the  future  before  computing  any  averages.  Details  of 
our  algorithm  will  be  presented  in  [ 1 6] .  Our  results  (see  Tables  1  and  2)  show 
that  while  the  variational  approximation  is  poor,  the  approximate  formula  (Eq.  7) 
is  quite  good.  Furthermore,  excellent  agreement  with  exact  results  is  obtained  by 
sampling  from  the  pure  |  o  ]  2  distribution. 
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Energy  Derivatives.  While  conventional  ab  initio  approaches  regularly 
compute  the  analytic  gradient  of  the  energy  with  respect  to  nuclear  coordinates 
in  order  to  determine  forces,  and  thereby  equilibrium  geometries  [17]  and  (by 
finite  difference  or  higher  analytic  derivatives)  harmonic  vibrational  frequencies 
[18],  only  finite-difference  approaches  have  been  implemented  in  QMC  [8].  In 
principle  there  is  no  reason  for  this  limitation.  To  compute  the  energy  derivative 
with  respect  to  a  nuclear  coordinate  p,  we  write 

d<E>f^__  d  I  jm)'j'T(E.)EL(R)dR 
dp  dp\  dR 


=  < 


dEL 

dp 


>, 


•  oo 


+  <if jEl  >,  - 

<p  dp  'oo 


<i|7>,  <El>, 

<p  dp  ;°o 


+  < 
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dp 


-E,  > ,  -< 


dV- 


dp 


>  <Er  > 


(8) 


The  second  equality  is  obtained  from  differentiation  using  the  chain  rule,  followed 
by  expression  of  the  resulting  ratios  as  averages  over  the  distribution  fx.  The 
derivative  dojdp  is  unknown:  it  is  however  possible  to  sample  it.  The  other 
terms  in  Eq.  8  may  be  evaluated  straight-forwardly  during  the  QMC  simulation. 
Rather  than  sampling  do/ dp,  as  a  first  approximation  we  may  take 
O  1  do  /  dp='i>  dty  T  /  dp.  This  turns  out  to  be  a  good  approximation  even  when 


^7-  is  only  of  moderate  accuracy  (e.g.  double- zeta  Hartree-Fock). 

Using  this  approach,  we  have  performed  calculations  on  H2  at  a  few  nuclear 
separations  [19].  Our  results  are  presented  in  Table  3,  and  Figures  1  and  2, 
where  they  are  compared  with  the  essentially  exact  work  of  Kolos  and  Wol- 
niewicz  [20],  as  well  as  with  the  results  of  conventional  ab  initio  approaches.  As 
can  be  seen,  QMC  is  competitive  with  Cl,  and  far  superior  to  Hartree-Fock. 
Attempting  to  construct  the  potential-energy  curve  of  H2  from  just  the  four 
QMC  energy  data  points  (all  of  which  are  exact  to  within  the  standard  error) 
leads  to  the  curves  of  Figure  1.  However,  the  additional  information  present  in 
the  first  derivatives  leads  to  a  Monte  Carlo  potential-energy  curve  (cf.  Fig.  2) 
indistinguishable  from  the  exact  one  [20]. 

Excited  States.  As  mentioned  in  Sect.  I,  work  thus  far  with  QMC  has 
been  limited  to  ground-state  potential-energy  surfaces  and  lowest-energy  states  of 
a  particular  symmetry  [4-7,14].  For  example,  we  have  calculated  the  energy  of 
the  first  excited  state  of  methylene  [7]  in  order  to  obtain  the  (until  recently) 
elusive  singlet-triplet  splitting.  This  was  the  first  molecular  QMC  calculation  of 
an  excited  state.  Our  results  there  were  in  excellent  agreement  with  the  most 
recent  experiments.  The  restriction  on  lowest  energy  states  of  a  symmetry  comes 
from  an  essential  feature  of  the  mapping  of  the  Schrodinger  equation  into  its 
diffusion  analog— namely,  that  time  in  these  two  equations  differs  by  a  factor  of  t  . 
This  means  that  the  expansion  of  a  time-dependent  molecular  state  vector  in 
energy  eigenfunctions  multiplied  by  exp (-iEt  /ft),  results  in  a  series  in  which  only 
the  lowest  energy  term  (i.e.  <p )  survives  at  large  t .  Thus  one  obtains  exponential 
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convergence  to  the  lowest  energy  eigenstate. 

If  'I' T  is  orthogonal  to  the  exact  lowest-energy  state,  one  can  see  from  Eq.  6 
that  convergence  will  be  to  the  next-lowe'-t  energy.  (Initially  4>  contains  a  super¬ 
position  of  states.)  This  fact  actually  is  used  even  in  the  calculation  of  molecular 
ground-state  energies,  as  Fermi  states  are  excited  states  (with  respect  to  the 
Boson  ground  state)  of  the  Schrodinger  equation.  By  choosing  T  to  be 
antisymmetric  with  respect  to  particle  exchange,  one  can  project  out  all  sym¬ 
metric  states.  A  similar  result  holds  for  calculations  of  different  symmetry  states 
of  a  given  molecule. 

When  studying  states  of  the  same  symmetry,  it  is  generally  not  possible  to 
find  a  trial  wave  function  exactly  orthogonal  to  all  the  lower-energy  states  of  that 
symmetry.  This  implies  that  convergence  will  ultimately  be  to  the  lowest-energy 
state.  However,  the  fixed- node  approximation  used  to  treat  the  Fermi  problem  is 
of  assistance  here  too.  In  the  fixed-node  approximation,  the  nodes  of  V  T  are 
used  to  divide  R  -space  into  separate  volume  elements.  The  Schrodinger  equation 
is  solved  separately  in  each  of  these  elements.  This  results  in  a  solution  of  the 
Schrodinger  equation  with  added  boundary  conditions.  Viewed  this  way,  the 
Fermi  problem  is  handled  by  forcing  the  generation  of  an  antisymmetric  state 
above  the  Bose  ground  state  through  the  placement  of  nodes  in  the  solution  o. 
In  like  manner,  other  excited  states  can  be  treated  approximately  by  imposing 
additional  nodes.  The  accuracy  of  the  approximation  will  depend  on  how  well 
these  nodes  are  placed. 


In  treating  excited  states,  traditional  ab  initio  methods  generate  wave  func¬ 
tions  with  the  correct  number  and  dimensionality  of  nodes.  Thus  such  wave 
functions  are  a  good  place  to  begin  in  searching  for  an  excited-state  ^  T .  We 
report  preliminary  results  achieved  with  QMC  in  this  way.  A  more  detailed 
report  is  in  preparation  [21].  Tables  4  and  5  report  our  results  on  respectively 
the  first  two  excited  state  of  the  He  atom  (which  are  of  the  same  symmetry  as 
the  ground  state)  and  of  the  first  excited  state  of  Hg  (which  is  of  a  different  sym¬ 
metry  than  the  ground  state).  In  the  case  of  He,  we  have  obtained  64%  of  the 
correlation  energy  for  the  ls2s  state  and  90%  of  the  correlation  energy  of  the 
ls3s  state.  Though  the  former  appears  low  (percentage-wise),  we  note  that 
our  total  energy  is  within  0.66  kcal/mol  of  the  experimental  energy,  which  is  gen¬ 
erally  regarded  as  chemical  accuracy.  For  we  note  that  there  is  some  basis-set 
dependence— at  least  with  such  small  bases.  Nevertheless,  a  fairly  simple  basis  set 
(double-zeta  plus  polarization)  yields  approximately  80%  of  the  correlation 
energy.  We  expect  that  better  results  will  be  obtained  through  the  use  of  better 
optimized  trial  functions. 

Improved  LiH.  In  test  runs  on  our  properties  and  derivatives  programs,  as 
well  as  in  tests  on  our  new  correlation  functions  [see  (1)  in  Appendix],  we  have 
been  using  the  LiH  molecule  as  a  sample  system.  In  the  course  of  these  runs,  we 
have  calculated  new  bounds  on  the  ground-state  energy  of  LiH.  In  fact,  our  new 
correlation  function  has  allowed  us  to  rapidly  reduce  our  variance  for  this  energy. 
Our  previous  LiH  results  [4]  were  better  than  any  ab  initio  calculations  available 
at  the  time.  Recently,  workers  using  Cl  techniques  have  put  a  large  effort  into 


showing  that  they  could  do  better  than  QMC.”  They  achieved  a  newr  energy 
bound,  lower  than  our  old  QMC  results.  However,  as  a  result  of  our  test  runs 
alone,  we  have  achieved  an  energy  which  is  essentially  the  exact  result  (Table  6). 
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Appendix 

Related  Work  Supported  by  DOE 

Trial  Wave  Functions.  For  use  in  QMC  one  wants  a  trial  function  which 
is  as  simple  as  possible,  since  it  will  require  repeated  evaluation  at  each  step  of 
the  random  walk.  Yet  one  wants  a  function  which  provides  accurate  results.  In 
principle,  since  QMC  solves  the  Schrodinger  equation,  one  should  obtain  accurate 
results  regardless  of  the  choice  of  'I'j-.  However,  as  we  noted  already,  for  Fermi 
systems  inaccurate  nodes  in  will  lead  to  a  small  error  when  the  fixed-node 
approximation  is  used.  Furthermore,  the  statistical  “noise’'  will  be  large  for  a 
poor  choice  of  ^7- . 

We  have  found  in  our  work  that  a  single  determinant  f  with  only  a 
double-zeta  basis  set  places  the  nodes  extremely  well  in  ground-state  calculations, 
as  determined  by  the  quality  of  the  computed  total  energies.  Increasing  the  basis 
set  beyond  double  zeta  appears  to  offer  insignificant  gain  in  either  accuracy  (i.e. 
the  fixed-node  error  does  not  noticeably  decrease)  or  precision  (the  statistical 
uncertainty,  for  equal  computing  time,  remains  essentially  unchanged).  In  prac¬ 
tice  we  have  included  an  electron-electron  Jastrow  factor  in  our  functions  'F  j  in 
order  to  reduce  statistical  fluctuations,  and  in  many  cases  we  have  also  included 
an  electron-nuclear  factor.  Neither  factor  affects  the  positioning  of  the  nodes, 
and  hence  the  fixed-node  energies  remain  unchanged.  However,  the  Jastrow 's  do 
have  an  important  effect  on  variance  reduction. 


By  optimization  of  a  small  number  of  parameters  (2  -  4),  it  appears  that  an 
order  of  magnitude  improvement  in  computer  time  can  be  achieved.  Addition¬ 
ally,  we  are  now  experimenting  with  other  forms  of  correlation  function  in  the 
hope  of  finding  further  variance  reduction.  In  contrast  to  the  usual  Jastrow 
forms  of 


^=exp{?T^-}’ 


(Ala) 


and 


(Alb) 

for  the  electron-electron  and  electron-nuclear  terms  respectively,  we  have  done 
some  experimentation  with  an  electron-electron  correlation  function  of  the  form 


na -«  exp(6r,j  +cr,j ))  .  (A2) 

O' 

A  similar  form  can  be  used  for  the  electron-nuclear  correlation  function.  As  with 
the  Jastrow  form,  only  the  variance  is  affected,  as  this  term  does  not  introduce 
any  nodes  of  its  own.  It  turns  out  that  this  form  is  considerably  more  flexible 
than  Eq.  (Al).  Although  an  extended  Jastrow  form  appears  to  offer  the  same 
flexibility,  it  appears  more  difficult  to  optimize.  Iterative  improvement  of  'I'p  is 
also  of  potential  benefit.  One  such  scheme— a  global  rescaling  of  tyT~ has  been 
described  by  us  [7].  A  more  powerful,  local  scheme  is  being  developed  currently. 


Electron  Affinity  of  Fluorine.  Although  properties  for  which  o  is  not  an 
eigenstate,  such  as  the  electronic  charge  distribution,  need  to  be  sampled  as 
described  in  point  (1)  above,  many  important  quantities  are  energy  related. 


These  quantities  can  be  computed  more  simply— though  two  different  energies 
must  be  calculated  very  accurately.  One  such  energy-related  quantity  is  the  elec¬ 
tron  affinity.  This  past  year  we  have  performed  an  accurate  QMC  calculation  of 
the  electron  affinity  of  fluorine  [22].  These  results— like  our  energies— exceed  the 
quality  of  the  best  variational  calculations,  and  are  in  excellent  accord  with  the 
recommended  experimental  value. 

Parallel  Computing.  Monte  Carlo,  as  alluded  to  earlier,  is  such  a  compu¬ 
tationally  intensive  activity  that  new  techniques  are  needed  if  one  wishes  to 
attack  large  systems  or  obtain  very  high  precision  (e.g.  better  than  99.99%).  One 
avenue  we  have  explored,  in  collaboration  with  the  Advanced  Computer  Archi¬ 
tecture  Laboratory  at  LBL,  is  the  use  of  parallel  computing  architectures  [23]. 
Briefly,  we  have  found  that  Monte  Carlo  can  be  readily  made  to  run  at  95% 
efficiency  on  an  8  processor  system.  We  explored  several  different  directions  for 
parallelizing  QMC,  as  well  as  load-balancing  techniques  to  keep  the  efficiency 
near  100%.  It  is  expected  that  a  slightly  restructured  parallel  code  can  run  at 
essentially  100%  efficiency  with  an  almost  unlimited  number  of  processors.  Thus, 
with  sufficient  memory,  precision  will  scale  as  the  square  root  of  the  number  of 
processors,  while  computing  time  for  a  fixed  precision  will  scale  inversely  almost 
linearly  with  processors. 
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Table  1.  Comparison  of  expectation  values  for  properties  of  H  .  The  properties 
studied  are  expectation  values  of  the  squared  distance  from  the  axis,  along  the 
H0  axis,  and  from  the  center  of^the  molecule  (in  bohr2).  The  electric  quadrupole 
moment,  Q  ,  (in  esu  cm“  X 10-"6)  can  be  derived  from  the  other  expectation 
values.  The  trial  function,  V  T ,  is  a  single- zeta-plus-bond  SCF  function,  multi¬ 
plied  by  electron-electron  and  electron-r.uclear  Jastrow  functions.  The  function. 
<t>,  is  the  exact  wave  function  in  this  case,  since  there  are  no  nodes  for  the 
ground-state  of  H  .  However,  a  small  bias  due  to  the  short-time  approximation 
may  be  present.  Here  the  time  step  is  r=0.01  au.  “Approximate  formula"  refers 
to  Eq.  7  of  the  text.  Statistical  uncertainty  in  the  last  significant  figures  is  shown 
in  parentheses. 


Method 

(<x2>  +  <y2>) 

<z2> 

1  9 

<r  > 

Q 

Best  Variational4 

1.554 

1.020 

2.574 

0.664 

<  *  T  |  A  |  *  T  > 

1.543(2) 

1.078(2) 

2.621(3) 

0.49(1) 

<^T  | A  |  o> 

1.534(5) 

! 

1.047(4) 

2.580(9) 

0.56(2) 

Approximate 

Formula 

■  -  '  "  1 

1.525(10) 

1.016(9) 

J 

2.539(18) 

0.63(5) 

<d  |  A  |  <5>> 

1.527(6) 

1.025(5) 

2.552(10) 

0.61(3) 

Exact15 

1.523 

1.023 

2.546 

0.61 

Table  2.  Comparison  of  values  obtained  for  the  electric  quadrupole  moment  of 
N0.  The  trial  function,  'I' T ,  is  a  double-zeta  SCF  function,  multiplied  by 
electron-electron  and  electron-nuclear  Jastrow  functions.  Time  steps  of  r=0.0025 
and  0.00125  au  are  used,  with  no  noticable  bias  present.  “Approximate  formula" 
refers  to  Eq.  7  of  the  text.  Statistical  uncertainties  are  indicated  in  parentheses. 


Method 


Q  (esucm  *X10  *6) 


Hartree-Fock‘ 


MCSCF 


<*7.  |Q  |  #r  > 


-2.19 


<*7-  |Q  \6> 


Approximate 

Formula 


Experiment 


Ref.  26. 


Table  3.  Comparison  of  energy  derivatives  of  H2  at  various  geometries  with 
standard  techniques.  The  QMC  trial  function  consists  of  a  double-zeta  SCF  wave 
function  multiplied  by  an  electron-electron  correlation  function  of  the  form 
J~[(l-a  exp(6ri; +cr,y )).  The  parameters  used  are  a  =0.48,  b  =0.5-1,  c  =0.33, 
O 

and  d=lA.  Time  steps  ranging  from  7=0.1  au  to  7=0.005  au  are  used.  The 
quoted  results  are  extrapolations  to  7=0.  Energies  are  in  hartrees  and  derivatives 
in  hartrees/bohr.  Statistical  uncertainties  are  indicated  in  parentheses. 


“  Ref.  28. 
b  Pulay  in  Ref.  17. 
c  Ref.  29. 
d  Present  work. 
f  Ref.  20. 


Table  4.  Comparison  of  the  energy  of  the  first  two  excited  *S  states  of  He  with 
SCF  and  experiment.  The  SCF  wave  function,  whose  energy  is  shown  in  the 
table,  is  used  as  the  QMC  'P  .  The  column  headed  tvcCE,  gives  the  percentage 
of  the  correlation  energy  recovered,  and  is  computed  relative  to  the  SCF  number 
in  the  first  row  of  the  table.  A  time  step  of  r=0.05  au  is  used.  Statistical  uncer¬ 
tainties  are  indicated  in  parentheses. 


ls2s 

ls3s 

Method 

E  j(M 

%CE 

„p) 

(kcal/mole) 

E2(h) 

%CE 

ME  ~E  exp) 

(kcal/mol) 

SCFa 

-2.14307 

0 

1.83 

-2.06036 

0 

0.5/  i 

QMC 

-2.14493(7) 

■ 

64 

0.66(4) 

-2.06119(7) 

90 

0.06(4) 

Experiment*1 

-2.14598 

IOC 

0.00 

-2.06128 

100 

000 
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Table  5.  Comparison  of  the  energy  of  the  first  excited  singlet  state  of  H2  ( B 
1E,f)  at  its  equilibrium  geometry  with  self-consistent  field  (SCF),  configuration 
interaction  (Cl)  and  exact  results.  For  QMC,  two  different  trial  functions  are 
used.  They  are  constructed  from  double-zeta  (DZ)  and  double-zeta-plus- 
polarization  (DZP)  SCF  functions.  The  column  headed  %  CE  gives  the  percen¬ 
tage  of  correlation  energy  recovered,  and  is  computed  relative  to  the  SCF  number 
in  the  first  row  of  the  table.  A  time  step  of  t-=0.01  au  is  used.  Statistical  uncer¬ 
tainties  are  indicated  in  parentheses. 
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-0.742  “ 

0 

9.4 
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41 
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QMC 

(DZP) 
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1.9 

Cl 
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0.9 

Exact 
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100 


0.0 
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Figure  Captions 

Figure  1.  Fit  of  the  potential  energy  curve  for  H„.  Two  different  fits  to  the 
four  energy  points  are  shown.  Neither  a  third-order  Legendre  polynomial  fit.  nor 
an  exponential  spline,  gives  an  adequate  representation  of  the  curve  compared 
with  the  exact  potential  of  Kolos  and  Wolniewicz. 


Figure  2.  Fit  of  the  potential  energy  curve  for  Hg  including  energy  derivatives. 
A  Hermite  polynomial  fit  to  the  data,  which  is  possible  with  the  additional  infor¬ 
mation  available  from  the  derivatives,  is  indistinguishable  from  the  exact  curve. 


Fiqure 


Fioure 


